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The exact nonlinear response of noninteracting (Bloch) electrons is examined within a nonequilib- 
rium formalism on the infinite-dimensional hypercubic lattice. We examine the effects of a spatially 
uniform, but time- varying electric field (ignoring magnetic-field effects). The electronic Green's 
functions, Wigner density of states, and time-varying current are aU determined and analyzed. We 
study both constant and pulsed electric fields, focusing on the transient response region. These 
noninteracting Green's functions are an important input into nonequilibrium dynamical mean field 
theory for the nonlinear response of strongly correlated electrons. 
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I. INTRODUCTION 

The linear- response theory of Kubc^ and Greenwood^i 
is an attractive approach to understand how electrons 
(in the solid state) interact with external electromag- 
netic fields. It can be used (in principle) to calculate 
general linear-response functions in systems that have ar- 
bitrarily strong electron correlations. Surprisingly, the 
linear-response regime for many bulk materials (espe- 
cially for parabolic band semiconductor aSif i .S i iS i i Z) and de- 
vices, holds for a wide range of electric field strengths. 

But there are a multitude of interesting nonlinear ef- 
fects in electric fields. Most electronic devices have a 
nonlinear current- voltage relation (transistors, Josephson 
junctions, etc.) and there is wide interest in nonlinear ef- 
fects in bulk materials as well (since it is the nonlinearity 
that often determines the ultimate performance). 

Devices are also becoming smaller and smaller. Semi- 
conductor processing line features are well below 100 nm, 
and there is significant research effort on nanoscale de- 
vices. In the latter potential difference of one volt 
produces an electric field on the order of 10^ V/cm for 
nanometer scaled devices. These fields are large enough 
for nonlinear effects to be important, if not critical, to de- 
termine the proper behavior in an external field. There 
also has been significant research performed on high en- 
ergy density pulsed laser experiments, where fields as 
high as 10^° V/cm can easily be attained over a short 
time scale. In that case, one drives the material out of 
equilibrium by the pulse, and studies how it relaxes back 
to an equilibrium distribution (as a means to determine 
relaxation times, etc.). 

There are few theoretical approaches to nonlinear ef- 
fects in solid-state systems. The formalism was developed 
independently by Kadanoff and BaymSi and Keldysh&iS 
in the early 1960s (Baymii and KeldyshiS have each writ- 
ten short historical accounts of their discoveries). These 
approaches include the effects of external fields to all or- 
ders and typically use perturbation theory to determine 
the effects of many-body interactionsi^. In the 1980s, 
Wilkins and collaborator3^i^i^iSiLiiiiSii& spent much ef- 
fort in developing these ideas further, and in examining 
nonlinear responses in finite dimensions. Here we extend 



that work to infinite-dimensional lattices, where we find 
many of the results for the electronic Green's functions 
can be determined analytically. Our formalism allows for 
an analysis of steady-state effects (like the Wannier-Stark 
laddersii) and of transient effects (like the response to 
a pulsed field). These noninteracting Green's functions 
are a necessary input to a complete nonlinear response 
dynamical mean field theory for strongly correlated elec- 
trons. We will present results for that work in a separate 
publication. 

The organization of this contribution is as follows: in 
Section II, we present the formalism for the nonlinear 
response, in Section III, we present our numerical results, 
and in Section IV, we present our conclusions. 



II. GREEN'S FUNCTIONS FOR BLOCH 
ELECTRONS IN AN EXTERNAL ELECTRIC 
FIELD 

The Hamiltonian for tight-binding electrons hopping 
on a hypercubic lattice (in the absence of any external 
fields) is 
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where ty is the Hermitian hopping matrix (chosen to 
bei^ tij = t* I2\fd for nearest neighbors as d ^ oo), and 
/I is the chemical potential. We shall consider the case 
when this system is coupled to an external electromag- 
netic field. An electromagnetic field is described by a 
scalar potential (/)(r, i) and a vector potential A(r, t) via 

E(r,t) = -V<^(r,t)-i^^ (2) 
c at 

for the electric field, with c the speed of light. We will 
use the Landau gauge where = to perform our cal- 
culations, so the electric field is described solely by the 
vector potential. This provides a significant simplifica- 
tion of the formalism for spatially uniform (but possibly 
time- varying) electric fields. 
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Unlike many time-dependent Hamiltonians, the effect 
of tfie vector potential is not easily described by adding 
a time-dependent piece to the Hamiltonian in addition 
to the time- independent piece in Eq. (Q. Instead, one 
uses the so-called Peierls' substitutioni^ for the hopping 
matrix: 



ti. exp 



A(r, t) ■ dr 



Ri 



(3) 



where R, is the spatial lattice vector associated with lat- 
tice site i (and similarly for site j) and e is the electric 
charge. Note that the Peierls' substitution is a simplified 
semiclassical treatment of the electromagnetic field (our 
vector potential is a classical, not quantum field) and we 
are ignoring dipole (and multipole) transitions between 
bands because we consider just a single-band model. In 
this case, the Hamiltonian of the noninteracting electrons 
coupled to an electromagnetic field becomes 



Hit) = ^ Uj exp 
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R, 



A(r,t) • dr 



The corresponding electric field becomes 



E(r,t) = - 



l9A(r,t) 
c dt ' 



i 

(4) 
(5) 



We will choose the vector potential in such a way that 
either the field is zero before t — and is then turned 
on, or the field becomes asymptotically small as i — > — oo 
and it is adiabatically switched on; in this way, the early 
time Hamiltonian is always given by Eq. (Q, and that 
will be used to establish the initial thermal equilibrium. 
The magnetic field has a complicated structure in infi- 
nite dimensions, because it involves the curl of the vec- 
tor potential, which would need to be defined correctly 
for the infinite-dimensional limit. Because we are in- 
terested in electric fields with weak spatial dependence, 
we shall assume the associated magnetic field is small 
enough that we can neglect it, even though we will allow 
the electric field to vary in time. This is an approxima- 
tion, because our electromagnetic fields no longer satisfy 
Maxwell's equations, unless the field is uniform in space 
and constant in time. This condition can be relaxed, per- 
haps by using a gradient expansion for the weak spatial 
dependence of the fieldsiS., but such an approach is cum- 
bersome in infinite dimensions. From now on, we neglect 
the spatial dependence of the vector potential (i.e., we 
are considering only spatially uniform but time-varying 
electric fields). 

It is convenient to introduce a momentum-space rep- 
resentation for the Hamiltonian, which becomes 



eA{t) . 

he 



(6) 



with Ck = '^3 exp[zRj-k] and cj^ = cj cxp[— iR^-k]. 
Note that the Hamiltonian in Eq. ||HJ) is a special time- 
dependent Hamiltonian, because it commutes with itself 



for all times \H{t) ,T-l{t')] = 0, which greatly simplifies 
the analysis of the time-dependent Green's functions de- 
veloped below. 

The expression for the time-ordered single-particle 
Green's function is defined to be 

g^(k,t,i') = -^(T(ck(t)4(t'))); (7) 

because of the special time dependence of the Hamil- 
tonian, this Green's function can be determined an- 
alytically. In Eq. 0, the operators are expressed 
in a Heisenberg picture, where the time dependence 
is 0{t) = exp[itn{t)]0 exp[-itn{t)] with H{t) deter- 
mined from Eq. I^J, the time ordering symbol T or- 
ders earlier times to the right (with a change of sign 
when two Fermionic operators are interchanged), and 
the angle brackets indicate a thermal averaging (O) = 
Tr[exp(-/37i:)0]/Tr[exp(-/3H)], with /3 = 1/T the in- 
verse temperature and the Hamiltonian being the field 
free (early-time) Hamiltonian from Eq. We directly 
solve for the Green's function by finding the time depen- 
dence of the momentum-dependent creation and annihi- 
lation operators, and then directly solve for the Green's 
function by taking the relevant expectation values and 
traces ^^ii^iSS. The starting point is to calculate the time 
dependence of the operators: 



eA(t) 
he 



)-m]4(0 



d , , z . eA(t) . T , , 

-rMi) = -T[e(k-^^)-A*]ck(t) 
dt h he 



which can be integrated to give 
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Ck(i) 



exp 
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— OO 
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sA(t) ^ 
he 

eA{t). 
he 



fi\dt 
— ii\dt 



(8) 
(9) 

[ (10) 
Ck.(ll) 



It is now easy to find the expression for the time- 
ordered Green's function by inserting the time depen- 
dence from Eqs. and into the definition of the 
Green's function in Eq. {Tj) to yield 



5^(k,t,0 



--e{t-t')cxp 

[l-/(6(k)-M)] 



sA(t) . 
he ' 



/i] dt 



+ -9(t'-t)exp 

h 

X /(e(k)-Ai), 



*r., eA(t). 

[e(k 7 — , 

he 



fi\dt 



(12) 



since the averages satisfy (cj^Ck) — /(e(k) — /i) and 
(ckc[) - [1 - /(e(k) - fi)] with fix) = + eMM] 
being the Fcrmi-Dirac distribution, and e(k) the band 
structure. 

In infinite-dimensional calculations, it is often impor- 
tant to also determine local properties, like the local 
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Green's function [gioc = X^k 5(^)1' '^^ local density of 
states (DOS). The Green's function in Eq. p2|l depends 
on both e(k) and e(k — eA/hc). Hence, the summation 
over momentum cannot be performed simply by intro- 
ducing an integral over the noninteracting DOS. Instead, 
the method of Mueller-Hartmann must be used^iiSSiSij to 
perform the integrations over the Brillouin zone and to 
extract the leading contributions as d ^ oo. The algebra 
is straightforward, but lengthy. The final result is 



X p(e) exp 



.e 1 

''hi 



X exp 
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he 



X e 
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ifi{t — t )/h 



he 



(13) 



where a denotes the component of the vector potential 
and p(e) = exp[—e'^/t*'^]/y/Trt*a'^ is the noninteracting 
DOS (and a is the lattice spacing). Note that in the limit 
A — > 0, this reduces to the well-known noninteracting 
Green's function on a hypercubic lattice. 

While the results of Eq. H13|l are completely gen- 
eral, they are quite cumbersome for calculations, and 
it is useful to consider some simpler limits. The eas- 
iest case to evaluate, which is what we consider for 
the remainder of this paper, is to examine the case 
where the vector potential lies along the (1,1,1,...) di- 
agonal [A{t) = ^(i)(l, 1, 1, ...)]. This choice simpli- 
fies the calculations significantly. In this case, the 
momentum-dependent Green's function in Eq. H12(l de- 
pends on just two macroscopic objects — the band struc- 
ture e(k) and an additional energy function e(k) = 
-t* linid^oo J2a sm{kaa) / Vd: 



exp 



i f\ eaAjt) 

- / {e(k)cos— 

h J^r he 



, . eaA(t) -. ,-" 
+ e(k) sin ' }dt 

he 



X ^--][e{t-t')-fie-^i)]. (w) 



Hence the local Green's function can be found by inte- 
grating over a joint density of statesSS 

P2(e, e) = 5] <5[e - eik)]S[e - e(k)l (15) 

k 

which yields 



gL{t,t') ^ de dep2{e,e)g'^{e,e,t,t'). 



(16) 



Using the techniques of Mueller-HartmanpSiiSSiS^ again, 
gives the following expression for the joint density of 
states: 



P2(e,e) 



1 



,i*2„.exp(-^-^). (17) 



Substituting the joint density of states of Eq. ((T7I) into 
Eq. (|16|l and integrating over e gives the final expression 
for the local Green's function: 



de[0{t~t')-f{e-fi)] 



, , . .e /■* - eaA{t) 

X p(e)exp[— i— / at cos — J 

h Jj, he 



X exp[— t* 



* . eaA{F)\ , 

dt sin — /4ft2 
t' he J 



X e 



ifi(t — t')/h 



(18) 



Of course, the result in Eq. H18() agrees with that of 
Eq. (|13|l when the vector potential lies along the diag- 
onal. 



III. NUMERICAL RESULTS 

We begin by studying the current density of the system 
in the presence of the electric field. The current operator 
is determined by the commutator of the polarization op- 
erator (defined by H = Ric|ci) with the Hamiltonian 
of the system. The expression for the a-component of 
the current-density operator has the following form: 



eat 
hVd 



^ sin I ka 



"-^) cW- (19) 



The expectation value of the ath component of the 
current can be easily calculated from the time-ordered 
Green's function in Eq. (|12|l in the limit t' t^: 



eat*'^ . / eaAa{t) 
Adnh \ he 



eaAa{t) 
he 



5^(k,i,t+), 



dfje - n) 
de — p(e), (20) 



where the summation over momentum is performed the 
same way as before. The total magnitude of the current 
density is just times this result, since each compo- 
nent along the diagonal is the same. In the limit of low 
temperature, we perform a Sommerfeld expansion, which 
gives 



(Jit)) 



eat*'^p{fi) . /'eaA{t)\ 
AVdnh \ he J 



(21) 



[with A{t) the value of the vector potential for each 
component]. Note that in the case of a constant field. 
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A{t) = —EctO{t) is a linear function of t, and the current 
is sinusoidal, even though the field is time-independent. 
This is the well-known Bloch oscillationf24 with a fre- 
quency LUBioch — eaE/h. Since we have no scattering, 
the system is a perfect conductor, but the periodicity 
of the lattice restricts the wavevector to lie in the first 
Brillouin zone which causes the oscillatory current. 

One can investigate a current-current correlation func- 
tion to determine a noise spectrum, but because the cur- 
rent is periodic, the noise profile would be just two delta 
functions for a constant field, and we won't learn any- 
thing interesting from examining the noise. 

It is interesting to note, that the current is nonzero for 
the case A{t) — const, which corresponds to the case of 
zero electric field. This is a consequence of the fact that 
the vector potential results in a shifting of the Fermi sur- 
face. In the case of an interacting system this current 
will be destroyed by interparticle scattering. In our case, 
a free-energy analysis will show that the lowest-energy 
state is the one without any current. There are a num- 
ber of analogies of the response of this system to the 
response of a superconductor (such as an ac response to 
a dc field, the presence of current-carrying states that do 
not disappear over time, etc.). All of these results are 
artifacts of the lack of scattering in the system. 

To find the resistivity of the system, we consider the 
case of a uniform static electric field (along the diagonal) 
of magnitude E^/d, which is turned on at i = 0, so that 
A{t) = -Ecte{t) [Aa{t) = ~Ecte{t)], and the potential 
along a path 6(1, 1, 1, ...)/Vd is equal toV — —EbVd (the 
length b is the distance over which we have a potential 
drop). The expression for the Ohm's law in the form V — 
jRa'^^^ (current density multiplied by the resistance-area 
product), gives the following expression for resistance- 
area product: 



Rd 



V 

3 



AirhEdb 



/ eaEt\ 

["ir) 



(22) 



The resistivity is defined to be 1/6 times the resistance- 
area product, in the linear-response limit of £^ ^ 0. 
Therefore, 



Pi 



in. resp. 



This result is proportional to d, as it should be because 
the conductivity is proportional to l/d in infinite dimen- 
sions. The correct resistivity is zero for a noninteracting 
system. Here we see that the linear-response resistivity 
in Eq. H23|) goes to zero in the limit of large time t — > oo. 

Let us estimate the linear response resistance of the 
ballistic metal from the expression in Eq. (|23|l . which 
can be finite because the linear-response resistance has a 
factor of b/t in it. For the ballistic metal the length b over 
which the electrons have moved in the time t should be 
b = vpt, with Wi? a suitable average of the Fermi velocity. 
This gives the resistance 

nun. resp. - ^2 ^d+lt*2 p^^^ ' ^^^J 



This expression corresponds to the Sharvin 
resistance2SiS& for a single-band model in infinite 
dimensions. In three dimensions, the Sharvin resistance 
is hjle^ divided by the number of channels, which is 
a Fermi surface factor multiplied by A.ti jk'^pArea. To 
compare with our formula, we must first note that we 
map the hopping integral onto the effective mass (for 
low electron filling) via 



h^^^d 



(25) 



and that af't* piy\i) = C is a constant of order one (pro- 
portional to {kFa)'^~^ for low filling). Therefore, 



4:TTmvFaVd h A\/d 

(X 



2e2 (kpay- 



(26) 



which has a Sharvin-like form (but appears to have the 
wrong dependence on kpa for d — 3; this most likely 
is an artifact of the problems with assuming a spherical 
Fermi surface in large dimensions, which is valid only for 
vanishing electron densities). 

We can also investigate the heat current carried when 
there is an electrical field present (but no temperature 
gradient), and we find that its average value vanishes at 
half filling, as expected, because the thermopower van- 
ishes at half filling, and we have no thermal gradients to 
directly drive a thermal current (in the general case, the 
energy part of the current vanishes, and the chemical po- 
tential piece will give a contribution of — /ij to the heat 
current). So heat transport is trivial unless one intro- 
duces a thermal gradient to the temperature, which we 
do not do here. 

Next we examine the spectral function and the density 
of states in the presence of a field. The time-dependent 
spectral function can be calculated from the retarded 
Green's function g^{t, t') = -{i/H)8{t-t'){{c{t),cHt')}) 
(with the operators expressed in a Heisenberg picture) 
using the Wigner coordinatesSl by introducing the aver- 
age time tave = (t+t')/2 and the relative time trei = t—t' 
variables. In this case, the spectral function as a func- 
tion of the average time (and Fourier transformed over 
the relative time) is equal to 



1_ 



yl(tai,e,k,Cj) = Ini / dtreie g^iKtavcUel) 



and the DOS is equal to 



1_ 



(27) 



^(W,W) = Im / dtreie'^*-^'gLitave,Uel)- (28) 



In general, the retarded Green's function can be found 
from the same technique used to calculate the time- 
ordered Green's function: first one introduces the time 
dependence of the Heisenberg operators, then one eval- 
uates the operator averages. Since the anticommutator 
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of two local creation and annihilation operators (or two 
operators in the momentum basis) is equal to one, we get 



, e(k) /■* eaA(t], 

X exp — — / rficos — 

n Jf., he 

e-(k) /•* . eoAit) 
'at sm — J 

he 



X exp[— I 



(29) 



for the momentum-dependent Green's function and 



, .e /■* - eaACt)^ 

X exp — I- / dicos — 

h Jf, he 



X exphi*2( j dtsin^^^^ /4ft'](30) 



for the local Green's function (using the t and t' coordi- 
nates). Note that these Green's functions have no tem- 
perature dependence, hence the spectral function and the 
DOS are independent of temperature. This is character- 
istic of a noninteracting system. 

The spectral function, in the absence of a field, is a 
delta function [v4(k,a;) = 6{uj — e(k) + //)]. When a field 
is turned on, the time dependence is no longer a pure 
exponential, so the spectral function deviates from the 
delta function, becoming a peaked function of nonvan- 
ishing width. In the limit where tave oo, the steady 
state is approached and the spectral function becomes 
a set of evenly spaced delta functions, since the Green's 
function becomes a periodic function in trei- 

The analysis for the local DOS is more complicated. 
Since the e dependence in Eq. is so simple, the inte- 
gral can be performed directly, with the result 



X exp 



j-*2 



4^2 



\i{j^ave •) trel ) | 



(31) 



where 



(it exp 



2A{t) 

he 



(32) 



In order to evaluate some numerical results, we first 
consider the case of a constant electric field turned on at 



t = 0. In this case, we get 

-^(iavej trel) ^( ^ave ^re//2)^( ^ai'e ^" ^rel /^^^rel 

/2) 

X [W + trel/2 + (1 - e^T^^*--*'-/^))-^] 
+ S{tave + trel/'2)S{~tave + trel/2) 

leah 



e{t 

ave 

h 

ieaE 



trel /2)9{tave — trel /2) (33) 



This result has some interesting properties. If i? ^ 0, 
then / — trel for all tave, and is a Gaussian in trei, 
which Fourier transforms to a Gaussian in frequency, i.e., 
it becomes the noninteracting DOS. There is an inter- 
esting scaling behavior. If we define tave = tave^-aE /h, 
trel = treieaE/h, and Co = cuh/eaE, then 



^itave 1 trel ) 



eaE 



I (tave 7 trel ) 5 



(34) 



with / a function independent of E. Hence 



9locitave: trel) 



':9{trei)e 



ifj.trci / eaE 



X exp 



t* 



Ae^a^E^ 



l-^ij'ave 1 ^rel ) | 



(35) 



and the DOS becomes 

A{tave,iC!) = -— Im 
TT 



dtrel e'"*"' gPoe [tave , Uel ) , (36) 



with the normalization chosen so J dujA{uj) = 1 (for eas- 
ier comparison of curves for different E). Hence we ex- 
pect the DOS to have the same shape as a function of 
uj (with a possible shift due to the chemical potential 
factor), but the amplitude of the oscillations grows as 
E increases [because of the minus sign in the exponent 
in Eq. H35() ]. But that turns out only to be true near 
Lu — 0. At other frequencies, the evolution with E is not 
always monotonic, because the DOS conserves total spec- 
tral weight, so there cannot be a monotonic evolution of 
the peaks at all frequencies. 

Note that the DOS satisfies two properties in equi- 
librium. The first is that the integral over frequency 
equals 1. The second is that the DOS is always posi- 
tive. The proof for the integral yielding 1 holds even 
in the nonequilibrium case, because the anticommutator 
of two Fermionic creation and annihilation operators at 
the same time is still one. The positivity does not hold, 
because the standard derivation, using the spectral rep- 
resentation, requires the Hamiltonian to be independent 
of time in order to be able to be used, and thereby prove 
the positivity. Indeed, the DOS in the presence of a field 
has regions where it is negative. 

It is interesting to consider the limit of large tave, i-e., 
tave OO, then we get the steady-state solution. We 
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take only the last term of I(tave, trei) in Eq. (|33|l because 
tave is always larger than t^eZ in this hmit. The Green's 
function becomes 



9loci^av 



oo, t. 



rel } 



:d{t 



rel } 



X exp 



t* 



, eaE 



(37) 



The Fourier transform of this is a set of delta functions, 
with different amplitudes, that are equally spaced in fre- 
quency, with a spacing eaE /h (since the Green's function 
is periodic in trei)- This is the famous Wannier-Stark 
ladder-'', expected for systems placed in an external elec- 
tric field. In the results plotted in Fig.^ the fact that the 
peaks at multiples of this frequency get larger, and grow 
in height as tave grows, indicates our results are showing 
the correct build-up to the steady state, but they will 
never get there until tave — > oo. It is no coincidence that 
this frequency is the same as the Bloch oscillation fre- 
quency. This discussion was first described in detail from 
the Green's function approach by Davies and Wilkinsi^. 
Note that the DOS is nonnegative in the steady state. 

We can calculate the weight of the delta functions by 
performing the Fourier series integral. The frequencies 
are NeaE/h, and the Fourier coefficient is 
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For our numerical results, we examine how the system 
approaches the steady state as the field is turned on. We 
work at half filling (/i = 0), where the DOS is symmetric; 
hence, we plot only the results for positive frequencies. 
The field needs to be large enough for our calculations 
to be able to see the nonlinear effects of the field on 
the DOS. For us, the numerical results can easily see 
effects on the DOS when eaE/h > 0.1. In Fig. ^ we 
plot results for eaE/h — 1. While it is true that the 
Green's functions for t and t' both less than zero are equal 
to their equilibrium (field-free) limit, the Wigner DOS 
feels the effect of the fields for all finite tave, because the 
integral over trei always includes some Green's functions 
with either t or t' larger than zero. We can see that 
significant "precursor" effects occur only for tave > — 2 
here, and the DOS develops significant oscillations before 
one can see the delta functions start to build up at the 
integer frequencies. 

We plot a close up of the region around = 1 in Fig.|21 
Note how a sharp peak develops as the average time in- 
creases, but there are significant oscillations near uj = 1 
whose amplitude decreases slowly as t^ve increases. 

In Fig. 13 we plot the DOS in the u) variable near ui — 



for tn 



100 and for five values of eaE/h (0.1, 0.3, 1.0, 



FIG. 1; Density of states A{tave,u;) for the noninteracting 
electrons with eaE/h = 1. Note how the DOS is essentially a 
Gaussian for tave < —2, but then develops large oscillations as 
tave increases. The DOS approaches a steady state for large 
time given by a set of delta functions, equally spaced by the 
Bloch oscillation frequency. The DOS is no longer positive 
once the field is turned on, but the integral does always equal 
1. 



3.0, and 10.0). This shows how the oscillations grow as 
E increases. For other integer values of uj, the evolution 
is not monotonic in the field strength E (for example, at 
u! — 1 the peak values increase with E for 0.1 < eaE/h < 
0.7 and then decrease for 0.7 < eaE/h < 10). 

In addition to the spectral function and the DOS, it 
is interesting to examine the distribution function. In 
equilibrium, the distribution function is a Fermi-Dirac 
distribution function, but the distribution function can 
change for nonequilibrium cases. In order to discuss dis- 
tribution functions, we need to define two more Green's 
functions — the so-called lesser and greater Green's func- 
tions. They are defined as g>{t,t') = -{i/h){c{t)c'^{t')) 
and g'^{t,t') = {i/h){c^{t')c{t)) (with the operators ex- 
pressed in a Heisenberg picture). These Green's func- 
tions can also be determined exactly for Bloch electrons, 
and their expressions are the same as those for the re- 
tarded Green's function in Eqs. and (|Sn|l . except the 
9{t — t') factor is replaced by — /(e(k) — ^) for g"^ and 
by [1 — /(e(k) — /i)] for . There are three cases for the 
distribution function that we can consider (the Wigner 
distribution, the quasiparticle distribution, and the local 
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0.9 0.95 1 1.05 1.1 
Frequency 

FIG. 2: Close up of the density of states A{tave, near u — 1 
for the noninteracting electrons with eaE/h — 1. Note how 
the DOS approaches a steady state for large time by devel- 
oping a sharp peak, but that there are significant oscillations 
near the sharp peak that decay slowly in time. 




-0.1 -0.05 0.05 0.1 
Reduced Frequency cj 



FIG. 3: Close up of the density of states A{iave,ijj) near 
Lj — for the noninteracting electrons with eaE/h = 0.1, 0.3, 
1.0, 3.0, and 10.0. Note how the peak in the DOS evolves as 
a function of the electric field. 



quasiparticle distribution) . The most often used distribu- 
tion function is the Wigner distribution function, defined 
to be 

The Wigner distribution function is always equal to the 



field-free Fermi-Dirac result fwigner{tave^^) = /(e(k) — 
^) for Bloch electrons. The quasiparticle distribution 
function is defined in analogy with the equilibrium result 
[g<(k, cj) = 27ri/(tj)A(k, Lj)] via 



f quasi {^ave ; k) 



1 Img<fta^e,k,a;) 

27r A{tave,K^^) 



(40) 



(note that the name quasiparticle distribution does not 
necessarily imply that there must be an underlying 
Fermi- liquid in the system). Since the only difference 
between the retarded Green's function and the lesser 
Green's function is the replacement of the theta func- 
tion by the Fermi-Dirac distribution (which does not 
depend on the time variables), the ratio of the two 
terms in Eq. H4U|) has an explicit factor of /(e(k) — fi). 
The Fourier transform of the numerator is over all trei, 
while the denominator is only over all positive trei- 
The integral I{tave,trei) is an odd function of trei [see 
Eq. (|32|) ]. which implies the numerator in Eq. (|40|l is 
27r/(e(k) — /i)A(to^e, k, w), and we find the quasiparticle 
distribution function is equal to the Fermi-Dirac distri- 
bution once again. The final distribution function to be 
defined is the local quasiparticle distribution function. 
This is 



loc ^, 1 Img<^(i 

quasi K^ave) 



locK'-avcT^ 



2n A{ta 



(41) 



This distribution function is nontrivial in a field, because 
the DOS and the lesser Green's function both have oscil- 
lations, but the zeros occur at different locations on the 
frequency axis, so the ratio in Eq. (|41|l can have signifi- 
cant oscillations. 

The calculation of the local quasiparticle distribution 
function is difficult because the presence of an /(e — /i) 
factor precludes us from performing the integral over e 
analytically; hence the numerical computations are more 
involved. We need to evaluate the integral 

9^{tave,trei) = / dep(e)/(e - ^) 



X exp 



i*2 



numerically. If eaE/h = 0, then this is just the Fourier 
transform of 2TTi f (cu) p{lu) , which gives the correct lesser 
function. If eaE/h ^ 0, then the Green's function has to 
be calculated numerically. Because the real part of the 
lesser Green's function is nonzero for a longer range in 
time than the imaginary part, the function will have 
more oscillations than the function. The results for 
a local quasiparticle distribution function are plotted in 
Fig. 21 As it follows from this figure, the local quasipar- 
ticle distribution function varies significantly from the 
equilibrium values as tave increases. This is because 
the 5^ Green's function has high frequency oscillations, 
which are not as strong in the DOS. The oscillations con- 
tinue as tave increases, but they become difficult to plot. 
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Of course the momentum-dependent quasiparticle distri- 
bution function is equal to the Fermi-Dirac distribution 
function for this problem. 




1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
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-5-4-3-2-1 12 3 4 5 
Frequency 

FIG. 4: Local quasiparticle distribution function fiodtave,!^) 
for noninteracting electrons with eaE/h = 1 and T — 0.1. 
Note how the local quasiparticle distribution function varies 
significantly from the equilibrium values as tave increases (the 
lowest panel is for tave ~ 2). This is because the Green's 
function has high frequency oscillations, which are not as 
strong in the DOS. The oscillations continue as tave increases, 
but they become difficult to plot. 

Finally, we study the time dependence of the DOS 
for the case of a sharp pulse during the period of time 
< t < ts- The second derivative of the vector poten- 
tial is proportional to the strength of the magnetic field 
(which we are neglecting), so we want to keep the sec- 
ond derivative small for the calculations to make sense. 
We choose the electric field to have the following time 
dependence: E{t) = E9{tE — t)9{t), which corresponds 
to a vector potential 

A{t) ^ -cEte{tE - t)e(t) - cEtEO{t - ts). (43) 

Note that these results are "singular" for the noninteract- 
ing case, because the final vector potential is a constant 
that can correspond to a current carrying state if the 
Fermi surface is shifted from the zone center. Because 
there is no scattering, such a current lives forever (but 
would decay in the presence of any scattering). Numeri- 
cal calculations show that the DOS deviates visibly from 



its equilibrium value during the times \t\ < treiax when 
the amplitude of the field is larger or on the order of t* ; 
the relaxation time treiax is on the order of the pulse time 

tE- 

The results of the calculations are presented in Fig. [S] 
for eaE/h — 1 (when eaE/h is much smaller than 1, the 
oscillations become hard to see). The nonequilibrium 
DOS shows oscillating behavior, which then decays as 
time increases. The results satisfy a symmetry relation, 
where the Wigner DOS is identical for tave and t'^yf, when 

We also consider the case of a smooth pulse with a 
smooth turn-on and turn-off of the electric field: A{t) = 
EctEexp{—t'^/t'^)/2 [which corresponds to an electric 
field E{t) = Et/tEexp{-t'^/t%)]. This field changes sign 
at i = and has it maximum amplitude at t = ±\/0.5. 
The Wigner DOS is symmetric in tave, so we only plot 
results for positive times in Fig. Note that at tave = 
the field has been on for a long time, so the result is far 
from a Gaussian. The amplitude of the peak in the DOS 
at (jj = is largest at tave — ±\/0.5, and decays rapidly 
for larger times. 

The proof of the symmetry relation for the Wigner 
DOS is rather straightforward to do. If the vector poten- 
tial A{t) has definite parity: A{—t) — zLA{t), then it is 

easy to see from Eq. that I{-tave,trel) = I{tave,trel) 

for even functions and I{—tave,trei) — I{tave,trei)* for 
odd functions. Since it is the modulus of / that enters 
into the calculation of A{tave:^jj): the DOS will satisfy the 
given symmetry rules. For the case of the constant-field 
pulse, we need to shift the time axis by t£;/2 and shift 
the vector potential by Mb/2 to have a vector potential 
that is odd in time. The shift of the vector potential has 
no effect on the modulus of /, since it contributes only a 
phase, while the shift in the time axis, is precisely what 
is needed to give the symmetry relation described above. 
For the Gaussian pulse, the vector potential is already an 
even function, and the symmetry relation follows directly. 

Note that we do not calculate the experimental probe 
of the reflectivity as a function of time after the initial 
pulse, because this system has no intrinsic scattering, so 
the optical conductivity is always a delta function peak 
at zero frequency, hence we would not learn anything 
interesting from such an exercise here. It would be inter- 
esting to probe such behavior in systems with intrinsic 
scattering mechanisms, to understand how the different 
relaxation mechanisms can be detected. 



IV. CONCLUSIONS 

We have studied the nonlinear response of Bloch elec- 
trons to an external time varying (but spatially homo- 
geneous) electric field by employing an exact nonequi- 
librium formalism on an infinite dimensional hypercubic 
lattice. We found that the current showed Bloch oscilla- 
tions, even when the electric field was constant in time, 
and we derived a form for the Sharvin-like resistance of 
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the system. 

The time-dependence of the DOS was calculated. We 
showed that it becomes a Wannier-Stark ladder for long 
times, but the transient evolution toward those discrete 
delta functions had a complex structure, that survives 
out to long times. We also examined a number of dif- 
ferent kinds of distribution functions, and showed that 
the most commonly chosen distribution functions re- 
tained the Fermi-Dirac form regardless of the strength of 
the electric field (but the local quasiparticle distribution 
shows complex oscillatory behavior). For pulsed fields, 
we saw the transient response build and then decay. The 
amplitude of the oscillations was proportional to the am- 
plitude of the electric field E for a wide range of field 
strengths, and we needed the field to be sufficiently large 
{eaE/h ^ t*) before they could be easily seen. Of course. 



the oscillations decay at times larger than the pulse time. 

These noninteracting Green's functions form the basis 
for a nonequilibrium dynamical mean field theory, which 
we are currently developing to study the nonlinear re- 
sponse of systems close to the Mott transition. Results 
of that work will appear elsewhere. 
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FIG. 5: Local DOS for the case of a sharp flat pulse with 
eaE/h = 1.0, Ie = 10.0, and various average times. The 
horizontal scale is the same in every panel, but the vertical 
scale changes in the different panels. By comparing figure (a) 
with figure (b), one can see that the response is identical for 

times tave and to„e that satisfy tave + t'ave = ts- 
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FIG. 6: Local DOS for the case of a smooth Gaussian pulse 
with eaE/h = 10.0, Ie = 1-0, and various average times. 
The results are completely symmetric between negative and 
positive average times, so we plot only the positive times here. 
Note how the oscillations are already strong at tave = 0, first 
increase slightly, then fade away as the average time increases. 



